import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance
np.random.seed(1)
p1 = np.random.multivariate_normal(mean=[0, 1], cov=[[1, 0], [0, 1]], size = 15)
p2 = np.random.multivariate_normal(mean=[2, 3], cov=[[1, 0], [0, 1]], size = 15)
x1 = p1[:,0]
y1 = p1[:,1]
x2 = p2[:,0]
y2 = p2[:,1]

# 均值
mean_val_1 = [np.mean(x1),np.mean(y1)]
mean_val_2 = [np.mean(x2),np.mean(y2)]
p = np.concatenate((p1, p2), axis=0)
sigma1=[np.var(x1),np.var(y1)]
sigma2=[np.var(x2),np.var(y2)]

# 离差
s1=sigma1*(len(x1)-1)
s2=sigma2*(len(x2)-1)

plt.scatter(p1[:, 0], p1[:, 1], color='r')
plt.scatter(p2[:, 0], p2[:, 1], color='b')


# VI = np.linalg.inv(np.cov([1,1],p1))
# distance.mahalanobis([1,1],p1,VI)


def mahal(point,x,y,mean_val):
    point=np.mat(point).reshape(2,1)
    mean_val=np.mat(mean_val).reshape(2,1)
    cov_mat=np.mat(np.cov(x,y))
    mT=np.transpose(point-mean_val)
    cov_i = np.linalg.inv(cov_mat)
    m=point-mean_val
    dis=mT*cov_i*m
    dis=np.sqrt(dis)
    # dis=np.sqrt(np.transpose(point-mean_val).reshape(1,2)*np.linalg.inv(cov_mat)*(point-mean_val))
    return dis

dis_mat=np.ones((300,300))

def calc_curv():
    for i in range(300):
        for j in range(300):
            dis1=mahal([i/100.0,j/100.0],x1,y1,mean_val_1)
            dis2=mahal([i/100.0,j/100.0],x2,y2,mean_val_2)
            dis_mat[i,j]=np.abs(dis1-dis2)

# 数据来源
# calc_curv()
# dis_mat.tofile('dis_mat.bin')
dis_mat=np.fromfile('dis_mat.bin').reshape(300,300)

am=dis_mat.argmin(axis=0)
am=np.array(am)
curv=[]
for i in range(len(am)):
    curv.append(i)

curv=np.array(curv)
print(curv)
plt.plot(am/100.0,curv/100.0)
plt.show()